rm(list=ls())
options(scipen=20)

setwd("") ##set your working directory.




library(foreign)
library(MASS)
library(arm)
library(lmtest)
library(car)
library(ggplot2)
library(stargazer)
library(effects)
library(dplyr)
library(rstatix)
library(lmtest)
library(sandwich)
library(psych) 




data<-read.csv(file.choose(), header=T) # Kim_Lim_JJPS_2022_Data.csv


names(data)
summary(data)


tapply(data$kof, data$group, mean) 


data$group_f<-as.factor(data$group)




#### effect of Korea video on kof ######


m1<-lm(kof~group_f+gen+age, data=data) # with basic control variables
summary(m1)
coeftest(m1, vcov = vcovHC(m1, type="HC1")) # robust standard error

m2<-lm(kof~group_f+gen+age+I(kov==1)+I(kofr==1), data=data) # with two additional control variables, you lose more observations.
summary(m2)

#### effect of Korea video on kop #####

m3<-polr(as.factor(kop) ~group_f+gen+age, data=data)

summary(m3, digits = 3)

m4<-polr(as.factor(kop) ~group_f+gen+age+I(kov==1)+I(kofr==1), data=data)

summary(m4, digits = 3)  



#### Now effect of China video on chf ######


m5<-lm(chf~group_f+gen+age, data=data)
summary(m5) ## group 3 insignificant  


m6<-lm(chf~group_f+gen+age+I(chv==1)+I(chfr==1), data=data)
summary(m6) ## group 3 insignificant  


#### finally effect of China video on chp ######

m7<-polr(as.factor(chp) ~group_f+gen+age, data=data)


summary(m7, digits = 3) ## group 3 insignificant


m8<-polr(as.factor(chp) ~group_f+gen+age+I(chv==1)+I(chfr==1), data=data)


summary(m8, digits = 3)  ## group 3 insignificant  





## Figures 


eff<-effect("group_f",m2, confidence.level = 0.90) ## 90% interval used


eff_m2<-as.data.frame(eff)


ggplot(eff_m2,aes(x=group_f, y=fit)) +
  geom_bar(position="dodge", stat="identity", color=c("gray","red", "gray"), fill=c("lightblue","lightblue", "lightblue"), alpha=0.7) +
  geom_errorbar(aes(ymin=lower, ymax=upper),
                size=.2,    
                width=.1,
                position=position_dodge(0.85)) + coord_cartesian(ylim = c(50, 75))+
  xlab("Treatment") +
  ylab("Feeling Thermometer (South Korea)") +
  scale_x_discrete( labels=c("Control Group", "Korean Aid", "Chinese Aid")) 




eff2<-effect("group_f",m4, confidence.level = 0.90)

prob<-as.data.frame(eff2)
prob$mean<-prob$prob.X4+prob$prob.X5 ## extract the probablity of respondents choosing somewhat or strongly agree (4 or 5)
prob$lower<-prob$L.prob.X4+prob$L.prob.X5
prob$upper<-prob$U.prob.X4+prob$U.prob.X5




ggplot(prob,aes(x=group_f, y=mean)) +
  geom_bar(position="dodge", stat="identity", color=c("gray","red", "gray"), fill=c("lightblue","lightblue", "lightblue") , alpha=0.7) +
  geom_errorbar(aes(ymin=lower, ymax=upper),
                size=.2,    
                width=.1,
                position=position_dodge(0.85))+  coord_cartesian(ylim = c(0.40, 0.75)) +  xlab("Treatment") + ylab("Intergovernmental Cooperation with South Korea\n(Somewhat/Strongy Agree)") + scale_x_discrete( labels=c("Control Group", "Korean Aid", "Chinese Aid")) 











##### China (Appendix) #### 



eff<-effect("group_f",m6, confidence.level = 0.90)
eff_m6<-as.data.frame(eff)


ggplot(eff_m6,aes(x=group_f, y=fit)) +
  geom_bar(position="dodge", stat="identity", color=c("gray","gray", "red"), fill=c("lightblue","lightblue", "lightblue"), alpha=0.7) +
  geom_errorbar(aes(ymin=lower, ymax=upper),
                size=.2,    
                width=.1,
                position=position_dodge(0.85)) + coord_cartesian(ylim = c(26, 36))+
  xlab("Treatment") +
  ylab("Feeling Thermometer (China)") +
  scale_x_discrete( labels=c("Control Group", "Korean Aid", "Chinese Aid")) 





eff2<-effect("group_f",m8, confidence.level = 0.90)

prob<-as.data.frame(eff2)
prob$mean<-prob$prob.X4+prob$prob.X5 ## extract the probablity of respondents choosing somewhat or strongly agree (4 or 5)
prob$lower<-prob$L.prob.X4+prob$L.prob.X5
prob$upper<-prob$U.prob.X4+prob$U.prob.X5



ggplot(prob,aes(x=group_f, y=mean)) +
  geom_bar(position="dodge", stat="identity", color=c("gray","gray", "red"), fill=c("lightblue","lightblue", "lightblue") , alpha=0.7) +
  geom_errorbar(aes(ymin=lower, ymax=upper),
                size=.2,    
                width=.1,
                position=position_dodge(0.85))+  coord_cartesian(ylim = c(0.10, 0.30)) +  xlab("Treatment") + ylab("Intergovernmental Cooperation with China\n(Somewhat/Strongy Agree)") + scale_x_discrete( labels=c("Control Group", "Korean Aid", "Chinese Aid")) 




### interaction models 




m2_int<-lm(kof~group_f*I(kov==1)+gen+age+I(kofr==1), data=data)
summary(m2_int) 

summary(m2_int, digits = 3)


m2_int2<-lm(kof~group_f*I(kofr==1)+gen+age+I(kov==1), data=data)
summary(m2_int2) 

summary(m2_int2, digits = 3)






data$val<-reverse.code(-1,data$val, mini=1, maxi=5)


m1_aov <-aov(val ~ group_f, data) # Appendix

m1_tukey<-TukeyHSD(m1_aov, which=c("group_f")) # Appendix





m2_int3<-lm(kof~group_f*val+gen+age+I(kov==1)+I(kofr==1), data=data)
summary(m2_int3) ## a rather weak interaction effect significant at 0.1 level


eff_int<-effect("group_f*val",m2_int3,xlevels=list(val=c(3,5)), confidence.level = 0.90)

eff_m2_int<-as.data.frame(eff_int)[c(1,2,4,5),]


p3<-ggplot(eff_m2_int,aes(x=as.factor(val), y=fit, fill=as.factor(group_f))) +
  geom_bar(position="dodge", stat="identity" , alpha=0.7) +
  geom_errorbar(aes(ymin=lower, ymax=upper),
                size=.2,    
                width=.1,
                position=position_dodge(0.85)) +  coord_cartesian(ylim = c(38, 85))+ 
  xlab("Perceived Value of Angkor as a World Heritage Site") +
  ylab("Feeling Thermometer (South Korea)") +
  scale_x_discrete( labels=c("Neutral", "Very Valuable"))+scale_fill_discrete(name = "Treatment", labels = c("Control", "Watched Video on Korean Aid"))










### Without any control variable (Appendix)

m9<-lm(kof~group_f, data=data)
summary(m9)
coeftest(m9, vcov = vcovHC(m9, type="HC1"))


Y<-data[data$kov>=1,]
summary(Y)
Y$group_f<-as.factor(Y$group)


m10<-lm(kof~group_f, data=Y)  #With only 451 observations (Reviewer comment)
summary(m10) 



m11<-polr(as.factor(kop) ~group_f, data=data)

summary(m11, digits = 3) 

m12<-polr(as.factor(kop) ~group_f, data=Y) #With only 451 observations (Reviewer comment)

summary(m12, digits = 3)  



m13<-lm(chf~group_f, data=data)
summary(m13) 


 
m14<-lm(chf~group_f, data=Y) #With only 451 observations (Reviewer comment)
summary(m14)



m15<-polr(as.factor(chp) ~group_f, data=data)


summary(m15, digits = 3) 


m16<-polr(as.factor(chp) ~group_f, data=Y) #With only 451 observations (Reviewer comment)


summary(m16, digits = 3)  




